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ABSTRACT 

One-step  methods    for  solving  integro-differential   equations 
are   studied  from  the  point   of  view  of  desiring  that  the  method  give 
good  accuracy  when  the  true   solution  asymptotically   goes  to   zero.      A 
test   equation   is   proposed  and  absolute   stability   is    investigated. 
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1.   INTRODUCTION 

Behavior  of  a  number  of  physical   and  biological   systems  may 
he   described  by   integro-differential   equations   of  the   form 

(1)  y    (x)    =   F    (x,    y(x),    J    K    (x,   t,    y(t))    dt ) . 

o 

A  familiar  example  is  the  population  equation  [2] 

2  ,X 

(2)  y  (x)  =  ay(x)  -  by  (x)  -  y(x)  J  K(x-t)  y(t)  dt , 


in  which  y  represents   a  population  and  x  represents  time.      The   integral 

factor  takes    into  account    "heredity",    or  in   general  the  past   history  of 

the   system.      Thus   in  applications   to  mechanics   such   integral   factors 

may  provide   for  component   fatigue. 

In  spite  of  their  importance,   only  very   recently  has  work 

been  done  on  the   development   of  efficient   and  accurate  numerical  methods 

for  solving  these  equations.      The   simplest   methods   proposed   [l,    3]   are 

the  one-step  methods,    analogous   to   Euler's   method  and  its    generalizations, 

in  which  the  left  hand  side  of   (l)    is    replaced  by   a  first-order 

difference  y ( x+h )   -  y ( x )      and  the   right   side   is    evaluated  at    x   (or, 
h 

more  generally,    is    some  convex  linear  combination  of  evaluations   at 
x  and  x+h).      Beginning  with   a  given   initial   value  y(o),   one  then   carries 
out   the   familiar  step-by-step  process   of  determining  an  approximate 
solution   at  points   nh,    n  =  1,    2,    ...    . 

Multistep  methods  have  also  been  proposed   [U,  6],   but   in  the 
interests   of  simplicity  we  have  restricted  ourselves  to   one-step  methods 
in  this  preliminary   study.      In  Section  2  we  look   at   some   stability 
questions    for  Euler's  method  and  in  Section   3  we   consider  possible 
advantages   of  using  a  modified  Euler's  method. 


2.      STABILITY  OF   EULER'S  METHOD 

In   a  recent  paper,    [l],    Cryer  and  Tavernini   discussed  an 
Euler' s -method  approach  to  the   computation  of  solutions   to   Volterra 
functional   differential   equations.      Under  very   general  hypotheses   they 
show  that   the  method   is   stable   in  the   sense  that  the  errors   are  bounded. 
Just   as    for  ordinary   differential   equations,    however,   one  would  expect 
that   a  stronger  stability  condition   is   required  in  order  to   compute 
valid  solutions   to   a  functional   differential   equation  with   a  decreasing 
exponential   solution.      We  here   investigate  this   question   in  a  particularly 
simple   context.  * 

Since   it   is    impossible  to   investigate  the  stability  properties 
of  a  method  as    it  would  apply  to   any   equation,    it   is    common   in  studying 
methods    for  ordinary   differential   equations   to   introduce  the  notion  of 
a   "test   equation" — a  simple   equation   for  which' the  method  is   easily 
analyzed  to  yield  information   on  how  the  method  will  work   for  a  large 
class   of  "similar"   equations.      The   standardly  used  test   equation   for 
ordinary  differential   equations    is 
(3)     y    (x)   +  ay  =  0        (a>o) 

with  initial   condition  y(o)   =  1«      Nothing  this   simple  will  provide   a 
valid  test    for  an   integro-differential   equation  method,    since   an  equation 
involving  both   an  integral   and  a  derivative   is   essentially   analogous   to 
a  second-order  differential   equation.      With  this   in  mind,   we   decided  on 
the  test   equation 

y      +   ay  +    (a2/U)      /     y(t)    dt   =   0 
(h) 

y  (o)  =  1. 


For  a>0,  this  has  the  solution 

(5)  y  =  exp  (-  ax/2)  [l  -  ax/2]. 

Although  this  is  not  the  simple  decreasing  exponential  obtained  for 

(3),  it  does  have  the  desired  testing  behavior  of  approaching  zero 

asymptotically  as  x-*  +  °°. 

To  solve  (h)   numerically,  we  define  x  =  nh  (n  =  0,  1,  ...) 

n 

and  y(x   )    =  y    .      Using  Euler's  method  together  with  the  trapezoidal 
rule   for  the  quadrature,   we  then  obtain  the   following  numerical 
approximation  to    (h) . 

y     =    (l  -   ah)y 


(6)  y2  =  6y1  -   (y/2)y 


0 


n-1 

yn+l  =   6yn  "  Y      I     yj    "    (Y/2)   y0  (n>l) 

where   6=l-ah-ah/8 

2  2 
and  y  =  a-  h   /h. 

In  order  to   study  error  propagation,  we  assume  that  yn   is   exact,    and 
that   an   error  E     is   introduced  into  y   .      This    error  then  propagates 
according  to 

(7)  E2  =   6E   , 

n-1 


E  ^    =   6E     -  y        I     E  . 
n+1  n  L       J 


Letting  A  =  ah,  one  defines  the  region  of  absolute  stability  of  the 

method  as  the  set  of  complex  values  of  A  for  which  E  •*■  0  as  n  ■*■  **». 

n 

Since  the  order  of  the  difference  equation  (T)  changes  with 

n,  we  use  an  indirect  approach  to  determine  the  behavior  of  E  .   Define 

n 

the  accumulated  error 
n 

T    =     Ye,. 

n  L     j 

Then,    from   (T),    T     satisfies   the   simple   3-term  recursion 

(8)  T  =    (1+6)    T     -    (y+6)    T  (n>l), 

n+I  n  n-1  — 

with  T     ■  0,    T     =  E . 


The   solution  to   equation    (5)    is 


(9)  T     =   El 


n 


[(l+5)2  -   My+6)]1/2 


{f,  -  4} 


where 


(10) 


((1+5)    ±    [(1+6)2   -   My+5)]1/2J 
=  1   -    (A/2)    -    (A2/l6)    ±    (A/l+)[A+A2/l6]l/2. 


q±   2   =   1/2     \  (1+6)    ±    [(1+6)"   -   U(y+6)] 


It    is   clear  that    if  both     q_      and    |q_      are  less  than  one,   then  T  ■*»  0   as 

1 nl '  ' ^2 '  n 

n-»-  °°  and  E     must    approach   zero  also.      If  we  now  write   E     =   T     -  T     ,    in 
n  n  n  n-1 

terms   of  q_    and  q     and  consider  the  various   possibilities,    it  becomes 

clear  that  the   conditions    I  q_  |<1,    I  q_  I  <1   are  also  necessary  for  E  ■+  0. 

1  TL  '  '    2 '  n 

Thus  the  region   of  absolute  stability   is   the   domain  in  which  both    |  q.  | 
and    |q-|    are  less  than  one.      Looking  at   real   values  of  A,  we  see  that 
the  condition   for  this    is   0<A<2.      Computer  trials    for   A   in  the  vicinity 
of  2   confirm  that    A   =  2   is   the   dividing  point  between  error  growth   and 


error  decay.      To   determine  the   situation   for  complex   A,   we  have  computed 
\n    |,     |q    I    on   a  grid  of  values    in  the  complex  plane.      The   region  on 
which  they   are  both  less   than  one   is   only   slightly   distorted  from  being 
the   disk    |  X  -  1 1    <  1    (See  the  Figure). 

For  comparison,    recall  that   the  region  of  absolute  stability 
for  Euler's  method  applied  to  the  ordinary  differential   equation    (3)    is 
precisely  the   disk    |A-1|<1,   where   again  A   =   ah.      Intuitively,    one  might 
feel   that   addition  of  an   integral   operator  term  and  the  accompanying 
summing  of  errors (see    (?))   would  have   a  serious   effect   on  the   region  of 
absolute   stability.      Our  analysis   shows  that  this    is  not  the   case. 

A  very   closely  related  question   is  whether  the  numerical 
solution  y      approaches   zero   as   n  ->  °°  as    it   ought.      This   question  may  be 
investigated  by   a  very   similar  procedure.      Define 
n 

(11)  \  =  Uy 

5=0 

Then,    using   (6)   with  y     =  1,   we  obtain 

(12)  Yn+1  =    (1+5  )Yn   -    (Y+<5)Yn_1  +    (y/2)  (n>l) 

with   initial   conditions   Y     =  1,   Y.    =   2  -   A.      This    is   an   inhomogeneous 

o       1 

difference  equation  with  solution 

(13)  Yn  "  2l^T  «3-2X-q2V  -  (3-2X-9l)q2n}  +  1/2. 

where  q  and  q  are  given  by  (10).   If  both  |q  |  and  | q? |  are  less  than 

one,  Y  -+■  1/2  as  n  ->  °°,  and  the  convergence  of  the  series 
n 

oo 

I     y.    implies  that  y .   ->  o   as   j   ■>  °°.      Thus  we   see  that,    as   commonly 
3=o     J  J 


Figure.   The  broken  line  denotes  the  "boundary  of  the  region  of 
absolute  stability  in  the  complex  X   plane.  For  comparison,  the  solid 
line  is  the  circle  \\  -   ll   =1. 


occurs    in  the  case  of  ordinary   differential  equations,    the  numerical 
solution  has   the  proper  asymptotic  behavior  within  the  region  of 
absolute  stability. 

In   summary,    it   appears   that   the   simple  Euler's  method  is 
perhaps   not   as   useful  as   a  universal  method  as   Cryer  and  Tavernini 
seem  to   suggest.      The  lack  of  absolute   stability   in  some   situations 
is   a  serious   problem.      Furthermore,   there   seems  to  be  no   straightforward 
way  to   determine   a  parameter   A  which,    for  a  general   integro- differential 
equation,    plays   roughly  the  same  role   as  the  ah   in  the  test   equation. 
That   is,    there   is  much  more  work  to  be   done  on  prediction  of  whether 
Euler's  method  is   likely  to  be  absolutely   stable   for  a  particular 
equation   and  a  particular  step  size  h. 


3.      GENERALIZED  EULER  METHODS 

In  this   section  we  briefly   consider  the  class   of 
one-step  methods  which  may  be  written 


(Ik)  yn+l 


-  =    (l-y)   F   ..    +  yF    , 
n+1  n 


where  0<y<  1  and  F     represents  the  right   side  of   (l)    evaluated  at   x  . 
n  n 

(When  F  involves  an  integral,   the   "evaluation"  involves   an  approximate 
quadrature,    of  course.)      For  y  =  1  we  have  the  simple  Euler's  method 
discussed  in  the  last  section.     When  y<l  the  method  becomes   implicit, 
but  the  increase  in  computational   complexity  may  be  offset  by  improved 
stability  properties    (i.e.    stability   for  larger  values   of  h).      In  the 
case  of  ordinary  differential  equations  this   is   certainly  true. 
Applying   (lM  to  the  test   equation    (3)   it   is   a  trivial   exercise  to   show 
that  the  method  is   "A-stable"    (absolutely  stable   for  all   X  such  that 
ReX>0)  whenever  0<y<— . 

It   seems   very  reasonable  that  this   result  will  also 
extend  to  integro-differential  equations.      In  an   attempt  to   do  this,  the 
approach  of  Section  2,  with  the  same  test   equation,   has  been  applied  to 
the  generalized  methods.      One  readily  obtains   a  stability  condition  of 
the   form    |q    |,    |q_|<l;   the   difficulty  arises   in  trying  to   determine   for 
what   values  of  X  these  conditions   hold.      Recall  that   in  the  last   section 
the  boundary  of  the  domain  of  absolute  stability  was    found  by  computing 
values  of    |q    |,    |qp|.      Such  a  computational  approach  has  no  hope  of 
demonstrating  absolute  stability  for  half  of  the  complex  plane.      Some 
other  method  must  be  devised.      When  X   is   real,  however,   it   is  possible 


to   show  analytically  that   absolute   stability  holds   for  all   A>0  when 
0<y<-. 

An  advantage  of  having  a  class   of  methods   depending  upon  a 
parameter    (y)    is  that   this   parameter  may   sometimes  "be  adjusted  so  as  to 
improve  the   numerical   solution  in  some   sense.      Thus   Liniger  and 
Willoughby    [5]  have  used  the   free  parameter  y   in  order  to   fit  the 
numerical   solution  to   exponentially-decreasing  behavior  and  thus   obtain 
improved  solutions  to   "stiff"  ordinary  differential  equations. 

We  have   done   some   experimentation  on  this    sort   of  idea  for 
the  case  of  integro-differential  equations.     A  difficulty  is  that  the 
solution    (5)   to  the  test   equation    (h)   is   initially  almost   a  straight 
line,   the   exponential   dominance   appearing  later.      Hence  it   is  hard  to 
see  how  a  y  might  be   chosen  to  mimic  the   step-by-step  behavior  of  the 
true   solution  throughout.      Believing  it   of  overriding  importance  to 
begin  with  a  good  numerical   solution,   we  have  opted  to  use   y  to   attempt 
to   get   an  exact   solution  at  the  point   x     =  h. 

Evaluating  the  exact   solution   (5)    at  h  we  get 

(15)  y(h)    =    (l-A/2)    exp(-A/2), 

while  the  numerical  method  (again  with  the  trapezoidal  rule  used  for 
the  quadrature)  yields 

f^\  -  (l-Ay-A2(l-y)/8) 

1  +  A(l-y)  +  A  (l-y)/8 

Equating  the  expressions  in  (15)  and  (l6),  we  find  that 

(IT)      (l-A/2)  exp(-A/2)  -  1  +  A 


A{(l-A/8)  -  (l-A/2)  (l+A/8)  exp(-A/2)} 


=  1-y 


10 


As    A   increases,    y  rapidly  decreases   to   zero    (when   X   is   about   k)    and 
turns   negative.      Hence  only  a  very  small  range  of  A  values  produce 
meaningful  values   of  the  parameter  y   from  equation    (IT).      There  seems, 
therefore,    to  be  very  little  value  in  pursuing  this   approach.      It   does 
appear  to   indicate,   however,    that   if  some  member  of  the  class  of 
generalized  Euler  formulas    (lU)   is  to  be  chosen  for  routine  use,   the 
best   choice  might  well  be  y  =  0    (backward  Euler),   both   for   initial 
accuracy   and  long-term  stability.      It   should  be  noted  that   all  members 
of  the   class   should  have  0(h    )   local   truncation  error  and  thus   should 
be   first   order.      (See    [l]    for  a  detailed  error  analysis   for  Euler's 
method    (y  =  l).)      It   might  well  be  worthwhile  to   develop  higher  order 
single-step  methods   for  integro-differential   equations. 


11 
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